
/*
 * BSD Licence:
 * Copyright (c) 2002 Herbert Sauro [ hsauro@cds.caltech.edu ]
 * Systems Biology Group [ www.sys-bio.org ]
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without 
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, 
 * this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright 
 * notice, this list of conditions and the following disclaimer in the 
 * documentation and/or other materials provided with the distribution.
 * 3. Neither the name of the <ORGANIZATION> nor the names of its contributors
 * may be used to endorse or promote products derived from this software
 * without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
 * DAMAGE.
 */

// Most of the original code for this class came from one of my Delphi classes
// which has the original comment:

// ***************************************************************************** 
// Basic Matrix Unit for Delphi, May 1996. Based on the original iMAP C matrix   
// library developed at Edinburgh, 1990.                                         
// Copyright c 1997-1997, Future Skill Software, Herbert M Sauro                 
// email: HSauro@fssc.demon.co.uk                                                
//
// Version 1.4, dated April 1997                                                 
// ***************************************************************************** 

// Other more minor influences for the C# version came from 
// the Java JAMA code, Paul Selormey's article at CodeProject  
// DotNetMatrix - Simple Matrix library for .NET
// http://www.codeproject.com/csharp/psdotnetmatrix.asp
// and finally Lutz Roeder's Mapack

  
using System;
using System.IO;
using System.Globalization;
using System.Collections;
using System.Runtime.Serialization;
using System.Runtime.Serialization.Formatters.Soap;

namespace Athena.Functions {
  
	/** <summary> Base Matrix Exception </summary> */
	public class EMatrixException : System.ApplicationException {
		/** <summary> </summary> */
		public EMatrixException (string message) : base (message) {
		}
	}

	/** <summary>
	 * Exception: Used to indicate that an operation was performed on incompatible matrices
	 * </summary> */
	class EMatrixSizeError : EMatrixException {
		public EMatrixSizeError (string str) : base (str) {
		}
	}

	/** <summary>
	 * Exception: Generated if the matrix was singular when the operation required a non-singular matrix
	 * </summary> */
	class ESingularMatrix : EMatrixException {
		public ESingularMatrix (string str) : base (str) {
		}
	}

	/** <summary>
	 * Exception: Matrix isn't square and it needs to be.
	 * </summary> */
	class ENonSquareMatrix : EMatrixException {
		public ENonSquareMatrix (string str) : base (str) {
		}
	}

	/** <summary>
	 * Exception: Indexing is out of bounds
	 * </summary> */
	class EMatrixIndexOutOfBounds : EMatrixException {
		public EMatrixIndexOutOfBounds (string str) : base (str) {
		}
	}

   
	/// <summary>
	/// This class defines a Matrix class for use in the .NET framework.
	/// Matrix stores data of type double. The matrix class represents the mathematical
	/// object 'matrix'. For consistency with other projects that might use this class,
	/// indexing for accessing matrix elements, starts at ZERO.
	/// 
	/// The class also supports named rows and columns. This is useful when rows and
	/// columns get swapped during numeric operations and a user needs to keep track
	/// of the rows and columns. Named rows and columns are particularly useful in some
	/// physical problems where the rows and columns correspond to physical entities.
	/// <code>
	/// Examples:
	///   Matrix A = Matrix (new double[,]{{1.2, 3.4}, {5.5, 6.7}));
	///   Matrix A = Matrix (3,4);
	///   Matrix A = Matrix.Random (10, 20);
	///   Matrix Z = Matrix.sub (X, Y);
	///   Matrix Z = Matrix.mult (X, y);
	///   X.transposeSelf();
	///   X.invertSelf();
	///   Y = X.invert();
	///   solution = A.solve (B)
	///   Y[1,1] = 1.2;
	///   Console.WriteLine (Y);
	///     
	/// </code>
	/// </summary>

	[Serializable()]
	public class Matrix : ICloneable {

		/// <summary> A matrix can be given a name</summary>
		public string name = "Matrix Name";  
		// Use jagged array because we can optimize 
		// certain operations such as row exchange
		// This has been made public so that it serializes in the XML
		/// <summary>
		/// Internal Jagged array representing the matrix
		/// </summary>
		public double[][] mxx;
  
		/// <summary> Array of column names </summary>
		public string[] columnNames;
		/// <summary> Array of row names </summary>
		public string[] rowNames;

		private static Random rnd = new Random(((int)DateTime.Now.Ticks));
		/// <summary>
		/// Static variable that holds the tolerance for deciding whether a number
		/// is zero or not. Used when calling isZero() and equals().
		/// </summary>
		public static double smallNumber = 1E-14;						

		// These are public so that they can be serialized
		public int r;
		public int c;
		
		// ---------------------------------------------------------------------
		/// <summary>
		/// Create an empty zero sized matrix. Could be resized later.
		/// </summary>
		// ---------------------------------------------------------------------
		public Matrix() : base() {
			r = 0;
			c = 0;
		}
		
		
		// Set up some default names for the rows and columns
		private void setDefaultLabels() {

			rowNames = new string[r];
			columnNames = new string[c];

			for (int i = 0; i < r; i++) 
				rowNames[i] = "R" + i.ToString();
			
			for (int i = 0; i<c; i++)
				columnNames[i] = "C" + i.ToString();

		}

		// ---------------------------------------------------------------------
		/// <summary>
		/// Creates an empty matrix of size 'r' by 'c' 
		/// 
		/// Usage:
		///   Matrix m = new Matrix(4,4);
		/// </summary>
		/// <param name="r"> Number of rows in the matrix</param>
		/// <param name="c"> Number of columns in the matrix</param>
		// ---------------------------------------------------------------------
		public Matrix(int r, int c) : this () {
			
			if ((r <=0) || (c <= 0))
				throw new EMatrixSizeError ("Arguments to Matrix constructor must be postive and greater than zero");

			this.r = r;
			this.c = c;
		   
			// Worth checking for negative values?
			this.mxx = new double[r][];
			for (int i = 0; i < r; i++)
				mxx[i] = new double[c];
			setDefaultLabels();
		}
		

		// ---------------------------------------------------------------------
		/// <summary>
		/// Copy constructor.
		/// Make a deep copy of the matrix specified in the argument.
		/// 
		/// Usage m = Matrix (myMatrix);
		/// </summary>
		/// <param name="cpy">Make a copy of this matrix</param>
		// ---------------------------------------------------------------------
		public Matrix(Matrix cpy) : this() {
			
			this.r = cpy.r;
			this.c = cpy.c;

			this.mxx = new double[cpy.r][];
			for (int i = 0; i < cpy.r; i++)
				mxx[i] = new double[cpy.c];

			for (int i=0; i<cpy.r; i++)
				for (int j=0; j<cpy.c; j++)
					this.mxx[i][j] = cpy.mxx[i][j];
			this.rowNames = new string[cpy.r];
			this.columnNames = new string[cpy.c];

			cpy.rowNames.CopyTo (this.rowNames, 0);
			cpy.columnNames.CopyTo (this.columnNames, 0);
		}


		// <summary>
		// GetHashCode
		// </summary>
		// <returns></returns>
		public override int GetHashCode() {
			return base.GetHashCode ();
		}

		// ---------------------------------------------------------------------
		/// <summary>
		/// Clone the current matrix. Clone performs a deep copy of the object.
		/// 
		/// Usage:
		///   A = B.Clone();
		/// </summary>
		/// <returns></returns>
		// ---------------------------------------------------------------------
		public virtual object Clone() {
			return new Matrix (this);
		}


		// ---------------------------------------------------------------------
		/// <summary>
		/// Copies a matrix object into a rectangular array
		/// 
		/// Usage:
		///	 double[,] A = m.ToArray();
		/// </summary>
		/// <returns>Returns a rectangular array, double],]</returns>
		// ---------------------------------------------------------------------
		public double[,] ToArray() {
			double[,] result = new double[this.r, this.c];
			for (int i=0; i<this.r; i++)
				for (int j=0; j<this.c; j++)
					result[i,j] = this[i,j];
			return result;
		}

		// ---------------------------------------------------------------------
		/// <summary>
		/// Create a matrix from a rectangular array argument 
		/// 
		/// Usage: Matrix m = Matrix (new double[,] { {1, 2, 3}, {4, 5, 6}, {7, 8, 9} });
		/// </summary>
		/// <param name="m">Rectangular Array</param>
		// ---------------------------------------------------------------------
		public Matrix (double[,] m) : this() {
			this.r = m.GetLength (0);
			this.c = m.GetLength (1);

			if ((r == 0) || (c == 0))
				throw new EMatrixSizeError ("Rectangular array to Matrix constructor must not be empty");

			this.mxx = new double[this.r][];
			for (int i = 0; i < this.r; i++)
				mxx[i] = new double[this.c];
			for (int i = 0; i < this.r; i++)
				for (int j = 0; j < this.c; j++)
					mxx[i][j] = m[i,j];
			this.setDefaultLabels();
		}


		// ---------------------------------------------------------------------
		/// <summary>
		/// Create an array from an array argument using the supplied row and 
		/// column labels. The size of the row and colume names arrays must match
		/// the size of the matrix.
		/// </summary>
		/// <param name="M">Rectangular Array</param>
		/// <param name="rowNames">List of row names</param>
		/// <param name="columnNames">List of column names</param>
		// ---------------------------------------------------------------------
		public Matrix (double[,] M, string[] rowNames, string[] columnNames) : this() {
			
			this.r = M.GetLength (0);
			this.c = M.GetLength (1);

			if ((r == 0) || (c == 0))
				throw new EMatrixSizeError ("Rectangular array to Matrix constructor must not be empty");

			if (rowNames.Length != this.r)
				throw new EMatrixSizeError ("The size of the row label array does not match the number of rows in the matrix");

			if (columnNames.Length != this.c)
				throw new EMatrixSizeError ("The size of the column label array does not match the number of columns in the matrix");

			this.mxx = new double[this.r][];
			for (int i = 0; i < this.r; i++)
				mxx[i] = new double[this.c];
			for (int i = 0; i < this.r; i++)
				for (int j = 0; j < this.c; j++)
					mxx[i][j] = M[i,j];

			rowNames = new string[this.r];
			columnNames = new string[this.c];
			rowNames.CopyTo (this.rowNames, 0);
			columnNames.CopyTo (this.columnNames, 0);
		}

		// ---------------------------------------------------------------------
		/// <summary>
		///  Create a matrix of size r by c filled with scalar value d 
		///  
		/// Usage:
		///   Matrix m = new Matrix (2, 2, Math.Pi);
		/// </summary>
		/// <param name="r">Number of rows</param>
		/// <param name="c">Number of columns</param>
		/// <param name="d">Scalar to fill the matrix</param>
		// ---------------------------------------------------------------------
		public Matrix(int r, int c, double d) : this() {
       
			if ((r <=0) || (c <= 0))
				throw new EMatrixSizeError ("Arguments to Matrix constructor must be postive and greater than zero");

			this.r = r;
			this.c = c;  
			this.mxx = new double[r][];
			for (int i = 0; i < r; i++)
				mxx[i] = new double[c];
			setDefaultLabels();

			for (int i=0; i<r; i++)
				for (int j=0; j<c; j++)
					this[i,j] = d;
		}

		
		/// <summary>
		/// Readonly property that returns the number of rows in the matrix
		/// </summary>
		public int rows {
			get { return r; }
		}
		
		/// <summary>
		/// Readonly property that returns the number of columns in the matrix
		/// </summary>
		public int cols {
			get { return c; }
		}


		// ---------------------------------------------------------------------
		/// <summary>
		/// Access Property
		/// <code>
		/// Usage: Matrix m = Matrix (2,3); 
		///        x = m[2,1]
		///        m[0,1] = 1.2;
		/// </code>code>
		/// </summary>
		// ---------------------------------------------------------------------
		public double this[int r, int c] {
			get { return mxx[r][c]; }
			set { mxx[r][c] = value; }
		}


		// Taken from Code Project:
		// Redhotglue C# ArrayObject
		// By Michael France.
		private string[] ReDimension(string[] OldArray, int arrNewLength) {
			// declare a larger array
			string[] NewArray = new string[arrNewLength];
			
			// determine if we are shrinking or enlarging
			int rowMax = 0;
			if (OldArray.Length < arrNewLength)
				rowMax = OldArray.Length;
			else
				rowMax = arrNewLength;
			
			// place values of old array into new array
			for(int row = 0; row < rowMax; row++) {
				NewArray[row] = OldArray[row];
			}
			
			return NewArray;
		}


		// Taken from Code Project:
		// Redhotglue C# ArrayObject
		// By Michael France.
		private double[,] ReDimension(double[,] OldArray,
			int arr1stDimLength,
			int arr2ndDimLength) {
			// declare a larger array
			double[,] NewArray = new double[arr1stDimLength, arr2ndDimLength];
			
			// determine if we are shrinking or enlarging
			const int FirstDimension = 0;
			const int SecondDimension = 1;
			int xMax = 0;
			int yMax = 0;
			// determine if we are shrinking or enlarging columns
			if (OldArray.GetUpperBound(FirstDimension) < (arr1stDimLength - 1))
				xMax = OldArray.GetUpperBound(FirstDimension) + 1;
			else
				xMax = arr1stDimLength;
			
			// determine if we are shrinking or enlarging rows
			if (OldArray.GetUpperBound(SecondDimension) < (arr2ndDimLength - 1))
				yMax = OldArray.GetUpperBound(SecondDimension) + 1;
			else
				yMax = arr2ndDimLength;
			
			// place values of old array into new array
			for(int x = 0; x < xMax; x++) {
				for(int y = 0; y < yMax; y++) {
					NewArray[x, y] = OldArray[x, y];
				}
			}		
			return NewArray;
		}


       	// ---------------------------------------------------------------------
		/// <summary>
		/// Resizes a matrix to a new dimension r,c. The resized matrix is returned,
		/// any data in the original matrix is preserved in the resized matrix. 
		/// The current matrix is not modified. Matrices can be resized to be smaller or 
		/// larger than the original.
		/// 
		/// Usage:
		///   Matrix m1 = Matrix.resize (m2, 6, 5);
		/// </summary>
		/// <param name="m">Matrix to resize</param>
		/// <param name="r">New number of rows</param>
		/// <param name="c">New number of columns</param>
		/// <returns>Returns the resized matrix</returns>
		// ---------------------------------------------------------------------
		public static Matrix resize (Matrix m, int r, int c) {
			
			Matrix cpy = new Matrix (m);
			cpy.resizeSelf (r, c);
			return cpy;
		}
		

		// ---------------------------------------------------------------------
		/// <summary>
		/// Nonstatic resize method. Resize the current matrix, preserving any data. 
		/// 
		/// Usage: 
		///   m.resizeSelf (5,8);
		/// </summary>
		/// <param name="r">New number of rows</param>
		/// <param name="c">New number of columns</param>
		// ---------------------------------------------------------------------
		public void resizeSelf (int r, int c) {
            
			// Make room for the longest rowNames array
			int longestArray = Math.Max (r, this.r);
			string[] tRow = new string[longestArray + 1];
			// Prepare the default names
			for (int i=0; i<longestArray + 1; i++) 
 				tRow[i] = "R" + i.ToString();
			// and copy over the current rowNames. Default names will be
			// left where the current rowNames were not copied to.
			for (int i=0; i<this.r; i++)
				tRow[i] = this.rowNames[i];

			// We do the same for the column names
			longestArray = Math.Max (c, this.c);
			string[] tCols = new string[longestArray + 1];
			// Prepare the default names
			for (int i=0; i<longestArray + 1; i++) 
				tCols[i] = "C" + i.ToString();
			// and copy over the current colNames. Default names will be
			// left where the current colNames were not copied to.
			for (int i=0; i<this.c; i++)
				tCols[i] = this.columnNames[i];

			// Resize the row and column labels
			this.rowNames = ReDimension (tRow, r); //this.rowNames = tRow;
			this.columnNames = ReDimension (tCols, c); //this.columnNames = tCols;

			// Convert the current matrix into a rectangular array
			double[,] oldArray = this.ToArray();

			// Resize the rectangular array
			double[,] newArray = ReDimension(oldArray, r, c); 
			// Convert the rectangular array into a matrix
			Matrix m = new Matrix (newArray);
		
			// ... and copy all the new information to the current matrix
			this.mxx = m.mxx;
			this.r = m.r;
			this.c = m.c;	
		}
		
        
		// ---------------------------------------------------------------------
		/// <summary>
		/// Make a deep copy of the current matrix. The destination matrix will
		/// hold the new copy and must have the appropriate dimensions.
		/// It is an error if the destination matrix is not the same size.
		/// 
		/// Usage: 
		///
		/// <code>
		/// src.copyTo (dest); 
		/// 
		/// dest.resize (src.r, src.c);
		/// src.copy (dest)
		/// </code>
		/// </summary>
		/// <param name="dest">Destination matrix</param>
		// ---------------------------------------------------------------------
		public void copyTo (Matrix dest) {
			if ((this.r != dest.r) || (this.c != dest.c))
				throw (new EMatrixException("source and destination matrices must be the same size"));

			dest.r = this.r;
			dest.c = this.c;

			dest.mxx = new double[this.r][];
			for (int i = 0; i < this.r; i++)
				dest.mxx[i] = new double[this.c];

			for (int i=0; i<dest.r; i++)
				for (int j=0; j<dest.c; j++)
					dest.mxx[i][j] = this.mxx[i][j];

			// Copy over column and row names
			this.rowNames.CopyTo (dest.rowNames, 0);
			this.columnNames.CopyTo (dest.columnNames, 0);
		}

		/// <summary>
		/// Save the matrix to a file as comma separated values.
		/// </summary>
		/// <param name="fileName">File name</param>
		public void saveCSV (string fileName) {
			TextWriter tw = new StreamWriter (fileName);
			try {
				for (int i=0; i<r; i++) {
					tw.Write ("{0}", this[i,0]);
					for (int j=1; j<c; j++)
						tw.Write (", {0}", this[i,j]);
					tw.WriteLine();
				}
			} finally {
				tw.Close();
			}
		}

		/// <summary>
		/// Save the matrix to a file with specified delimiter between values.
		/// </summary>
		/// <param name="fileName">File name</param>
		/// <param name="delimiter">Character to use to separate columns of data</param>
		public void saveCSV (string fileName, char delimiter) {
			TextWriter tw = new StreamWriter (fileName);
			try {
				for (int i=0; i<r; i++) {
					tw.Write ("{0}", this[i,0]);
					for (int j=1; j<c; j++)
						tw.Write (delimiter + " {0}", this[i,j]);
					tw.WriteLine();
				}
			} finally {
				tw.Close();
			}

		}

		/// <summary>
		/// Save a matrix as XML to a file.
		/// </summary>
		/// <param name="fileName">Name of file, include full path</param>
		public void saveXML (string fileName) {

			System.IO.FileStream fs = new System.IO.FileStream(fileName, System.IO.FileMode.Create) ;                        
			//System.Xml.Serialization.XmlSerializer xmlSer = new System.Xml.Serialization.XmlSerializer(this.GetType ());                    
			SoapFormatter xmlSer = new SoapFormatter();
			xmlSer.Serialize(fs, this);
			fs.Close();                                                          
		}


		/// <summary>
		/// Load a matrix from a file previously saved using SaveXML()
		/// </summary>
		/// <param name="fileName">Name of file, include full path</param>
		/// <returns>Returns the loaded matrix to the caller</returns>
		// Load a matrix from an XML file
		public static Matrix loadXML(string fileName) {
			//Open the XML file
			System.IO.FileStream fs = new System.IO.FileStream(fileName, System.IO.FileMode.Open) ;                         
			SoapFormatter xmlSer = new SoapFormatter();
			Matrix result = (Matrix) xmlSer.Deserialize (fs);
			fs.Close();
			return result;
		}


		// ---------------------------------------------------------------------
		/// <summary>
		/// Augment the current matrix with the matrix in the method argument
		/// <para>result = m1 | m2</para>
		/// m2 must be square and have the same number of rows as m1
		/// 
		/// Usage:
		///   m1.augment (m2)
		/// </summary>
		/// <param name="m">Argument to agument</param>
		// ---------------------------------------------------------------------
		public void augment (Matrix m) {
			if ((this.r != m.r) && (! m.isSquare()))
				throw (new EMatrixException("source and destination matrices must be compatible to augment"));

			int originalSize_c = this.c;
			this.resizeSelf (this.r, this.c + m.c);
			for (int i = 0; i < this.r; i++)
				for (int j = 0; j < m.c; j++)
					this[i, originalSize_c + j] = m[i,j];
			for (int i = 0; i < m.c; i++)
				this.columnNames[originalSize_c + i] = m.columnNames [i];
		}

		// ---------------------------------------------------------------------
		/// <summary>
		/// Checks if the matrix contains all zero elements. Values less than 1E-14 
		/// are assumed to zero. This value is set by the static variable, smallNumber.
		/// </summary>
		/// <param name="m">Matrix to test</param>
		/// <returns>True if matrix is a zero matrix</returns>
		// ---------------------------------------------------------------------
		public static bool isZero (Matrix m) {
			for (int i=0; i<m.r; i++)
				for (int j=0; j<m.c; j++)
					if (Math.Abs (m[i,j]) > Matrix.smallNumber)
						return false;
			return true;
		}


		// ---------------------------------------------------------------------
		/// <summary>
		/// Returns true if matrices x and y are the same size. 
		/// 
		/// Usage:
		/// if (sameDimensions (m1, m2)) { ... }
		/// </summary>
		/// <param name="x"></param>
		/// <param name="y"></param>
		/// <returns></returns>
		// ---------------------------------------------------------------------
		public static Boolean sameDimensions (Matrix x, Matrix y) {
			if ((x.r == y.r) && (x.c == y.c))
				return true;
			else return false;
		}


		// ---------------------------------------------------------------------
		/// <summary>
		/// Returns true if the matrix is square  
		/// 
		/// Usage:
		///    if (m.isSquare()) { ... }
		/// </summary>
		/// <returns></returns>
		// ---------------------------------------------------------------------
		public Boolean isSquare () {
			if (r == c)
				return true;
			else return false;
		}

		// ---------------------------------------------------------------------
		/// <summary>
		/// Static method for generating an Identity matrix of a given size. 
		/// 
		///	Usage:
		///   m = Matrix.Identity(5);
		/// </summary>
		/// <param name="n">Size of identity matrix</param>
		/// <returns>Returns the identity matrix</returns>
		// ---------------------------------------------------------------------

		public static Matrix identity(int n) {
			Matrix m = new Matrix (n,n);
			for (int i=0; i<n; i++)
				m[i,i] = 1.0;
			m.setDefaultLabels();
			return m;
		}

		// ---------------------------------------------------------------------
		/// <summary>
		/// Zero the current matrix, all elements are set to 0.0.
		/// </summary>
		// ---------------------------------------------------------------------
		public void zero() {

			for (int i=0; i<r; i++)
				for (int j=0; j<c; j++)
					this[i,j] = 0.0;
		}

		// ---------------------------------------------------------------------
		/// <summary>
		/// Static method to generate a diagonal matrix with a given scalar value.
		/// 		
		/// Usage:
		///   m = Matrix.diagonal(4, 1.0);
		/// </summary>
		/// <param name="n">Size of matrix</param>
		/// <param name="value">Value in diagonal</param>
		// ---------------------------------------------------------------------
		public static Matrix diagonal (int n, double value) {
			Matrix m = new Matrix (n,n);
			for (int i=0; i<n; i++)
				m[i,i] = value;
			m.setDefaultLabels();
			return m;

		}

		// ---------------------------------------------------------------------
		/// <summary>
		/// Static method to create a matrix filled with uniformly generated random numbers.
		/// 
		/// Usage: 
		///   m = Matrix.RandomMatrix (4,5)
		/// </summary>
		/// <param name="r">Number of rows in the matrix</param>
		/// <param name="c">Nunber of columns in the matrix</param>
		/// <returns>Returns a matrix filled with uniformly distributed random numbers</returns>
		// --------------------------------------------------------------------- 
		public static Matrix Random (int r, int c) {
			Matrix m = new Matrix (r, c);
			m.name = "rm";
			for (int i=0; i<r; i++)
				for (int j=0; j<c; j++)
					m[i,j] = rnd.NextDouble();
			m.setDefaultLabels();
			return m;
		}
		
			
     	// ---------------------------------------------------------------------
		/// <summary>
		/// Write the matrix contents to a string. 
		/// 
		/// Usage:
		///   Console.WriteLine (m.ToString());
		/// </summary>
		// ---------------------------------------------------------------------

		public override string ToString() {
			
			using (StringWriter writer = new StringWriter(CultureInfo.InvariantCulture)) {

				writer.Write ("     ");
				for (int i=0; i<c; i++)
					writer.Write ("{0,-8}", this.columnNames[i]);
					writer.WriteLine();

				for (int i=0; i<r; i++) {
					writer.Write ("{0,-4}", this.rowNames[i]);
					for (int j=0; j<c; j++)
						writer.Write ("{0,7:F4} ", this[i,j]);
					writer.WriteLine();
				}
			return writer.ToString();
			}
		}
	    

		// ---------------------------------------------------------------------
		/// <summary>
		/// Swap two rows in the current matrix.
		/// </summary>
		/// <param name="r1">First row</param>
		/// <param name="r2">Second row</param>
		// ---------------------------------------------------------------------
		public void swapRows (int r1, int r2) {
			
			if (((r1 < 0) || (r1 > this.r)) || ((r2 < 0) || (r2 > this.r)))
			   throw new EMatrixIndexOutOfBounds ("Matrix index out of bounds in swapRows");

			double[] tr = new double[this.c];

			tr = this.mxx[r1];
			this.mxx[r1] = this.mxx[r2];
			this.mxx[r2] = tr;

			string tmp = (string) this.rowNames[r1];
			this.rowNames[r1] = this.rowNames[r2];
			this.rowNames[r2] = tmp;
		}

   
		// ---------------------------------------------------------------------
		/// <summary>
		/// Swap two columns in the current matrix.
		/// </summary>
		/// <param name="c1">First column</param>
		/// <param name="c2">Second column</param>
		// ---------------------------------------------------------------------
		public void swapColumns (int c1, int c2) {

			Matrix tr = Matrix.transpose (this);
			tr.swapRows (c1, c2);
			tr.transposeSelf ();
			tr.copyTo (this);
		}


		/// <summary>
		/// Delete a specified row from the current matrix.
		/// </summary>
		/// <param name="r1">Row to delete</param>
		public void deleteRow (int r1) {

			if ((r1 < 0) || (r1 >= r))
				throw new EMatrixIndexOutOfBounds ("Row index out of bounds in deleteRow");

			int rowCount = 0;
			Matrix result = new Matrix (this.r - 1, this.c);
			for (int i=0; i<this.r; i++)
				if (i != r1) {
					result.rowNames[rowCount] = this.rowNames[i];
					for (int j=0; j<this.c; j++)
						result[rowCount,j] = this[i,j];
					rowCount++;
				}
			for (int i=0; i<this.c; i++)
				result.columnNames[i] = this.columnNames[i];
			this.resizeSelf (result.r, result.c);
			result.copyTo (this);
		}

		/// <summary>
		/// Insert a zeroed row into the current matrix before the row 
		/// specified in the argument.
		/// </summary>
		/// <param name="r1">Insert an empty row BEFORE this row</param>
		public void insertRow(int r1) {
			
			if ((r1 < 0) || (r1 > this.r))
				throw new EMatrixIndexOutOfBounds ("Row index out of bounds in insertRow");

			// Prepare space for the new row, including a new row label
			double[][] m = new double[this.r + 1][];
			for (int i=0; i<this.r + 1; i++)
				m[i] = new double[this.c];

			string[] tRow = new string[this.r + 1];

			// Copy over top section
			for (int i=0; i<r1; i++) {
				m[i] = this.mxx[i];
				tRow[i] = this.rowNames[i];
			}
			// Copy over the new row
			m[r1] = new double[this.c]; // empty row
			tRow[r1] = "R" + (this.r+1).ToString();

			for (int i=r1; i<this.r; i++) {
				m[i+1] = this.mxx[i];
				tRow[i+1] = this.rowNames[i];
			}
			this.mxx = m;
			this.rowNames = tRow;
			this.r++;
		}

		/// <summary>
		/// Insert a zeroed row into the current matrix before the row 
		/// specified in the argument. The new row will be named using 
		/// the supplied name argument.
		/// </summary>
		/// <param name="r1">Insert an empty row BEFORE this row</param>
		/// <param name="name">New label to use in the row</param>
		public void insertRow(int r1, string name) {
			
			if ((r1 < 0) || (r1 > this.r))
				throw new EMatrixIndexOutOfBounds ("Row index out of bounds in insertRow");

			// Prepare space for the new row, including a new row label
			double[][] m = new double[this.r + 1][];
			for (int i=0; i<this.r + 1; i++)
				m[i] = new double[this.c];

			string[] tRow = new string[this.r + 1];

			// Copy over top section
			for (int i=0; i<r1; i++) {
				m[i] = this.mxx[i];
				tRow[i] = this.rowNames[i];
			}
			// Copy over the new row
			m[r1] = new double[this.c]; 
			tRow[r1] = name;

			// Copy over the bottom section
			for (int i=r1; i<this.r; i++) {
				m[i+1] = this.mxx[i];
				tRow[i+1] = this.rowNames[i];
			}
			// update the matrix with the new information
			this.mxx = m;
			this.rowNames = tRow;
			this.r++;
		}


		/// <summary>
		/// Insert a zeroed column into the current matrix before the column 
		/// specified in the argument.
		/// </summary>
		/// <param name="c1">Column to insert</param>
		public void insertColumn (int c1) {
			insertColumn (c1, "C" + (this.c+1).ToString());
		}


		/// <summary>
		/// Insert a zeroed column into the current matrix. The new
		/// column will be named with the supplied name argument.
		/// </summary>
		/// <param name="c1">Column to insert</param>
		/// <param name="name">Name to use to label column</param>
		public void insertColumn (int c1, string name) {

			if ((c1 < 0) || (c1 > this.c))
				throw new EMatrixIndexOutOfBounds ("Column index out of bounds in insertColumn");

			Matrix m = new Matrix (this);
			m.transposeSelf();
			m.insertRow (c1, name);
			m.transposeSelf();
			this.resizeSelf (m.r, m.c);
			m.copyTo (this);
		}


		/// <summary>
		/// Delete a specified column from the current matrix.
		/// </summary>
		/// <param name="c1">Column to delete</param>
		public void deleteColumn (int c1) {

			if ((c1 < 0) || (c1 >= c))
				throw new EMatrixIndexOutOfBounds ("Column index out of bounds in deleteColumn");

			Matrix m = new Matrix (this);
			m.transposeSelf();
			m.deleteRow (c1);
			m.transposeSelf();
			this.resizeSelf (m.r, m.c);
			m.copyTo (this);
		}



		/// <summary>Determines wheather two instances are equal.</summary>
		public override bool Equals(object obj) {
			return equals(this, (Matrix) obj);
		}

		// ---------------------------------------------------------------------
		/// <summary>
		/// Checks whether two matrices are equal. Two numbers are assumed to be
		/// different if their difference (a-b) is greater than the static variable
		/// smallNumber (default 1E-14).
		/// 
		/// Usage: 
		///   if (Matrix.equals (m1, m2)) ....
		/// </summary>
		/// <param name="m1">First Matrix</param>
		/// <param name="m2">Second Matrix</param>
		/// <returns>True if the matrices are equal</returns>
		// ---------------------------------------------------------------------
		public static bool equals (Matrix m1, Matrix m2) {
			
			if (sameDimensions (m1, m2)) {
				for (int i=0; i<m1.r; i++) {
					for (int j=0; j<m1.c; j++)
						if (Math.Abs (m1[i,j] - m2[i,j]) > smallNumber)
							return false;
				}
				return true;
			} else
				return false;
		}

		// ---------------------------------------------------------------------
		/// <summary>Returns the One Norm for the current matrix.</summary>
		/// <value>The maximum column sum.</value>
		// ---------------------------------------------------------------------
		public double Norm1 {
			get {
				double f = 0;
				for (int j = 0; j < c; j++) {
					double s = 0;
					for (int i = 0; i < r; i++) {
						s += Math.Abs(this[i,j]);
					}

					f = Math.Max(f, s);
				}
				return f;
			}		
		}

		// ---------------------------------------------------------------------
		/// <summary>Returns the Infinity Norm for the current matrix.</summary>
		/// <value>The maximum row sum.</value>
		// ---------------------------------------------------------------------
		public double NormInf {
			get {
				double f = 0;
				for (int i = 0; i < r; i++) {
					double s = 0;
					for (int j = 0; j < c; j++)
						s += Math.Abs(this[i,j]);
					f = Math.Max(f, s);
				}
				return f;
			}
		}

		private static double Hypotenuse(double a, double b) {
			if (Math.Abs(a) > Math.Abs(b)) {
				double r = b / a;
				return Math.Abs(a) * Math.Sqrt(1 + r * r);
			}

			if (b != 0) {
				double r = a / b;
				return Math.Abs(b) * Math.Sqrt(1 + r * r);
			}

			return 0.0;
		}

		// ---------------------------------------------------------------------
		/// <summary>Returns the Frobenius Norm for the matrix.</summary>
		/// <value>The square root of sum of squares of all elements.</value>
		// ---------------------------------------------------------------------
		public double FrobeniusNorm {
			get {
				double f = 0;
				for (int i = 0; i < r; i++) {
					for (int j = 0; j < c; j++) {
						f = Hypotenuse(f, this[i,j]);
					}
				}

				return f;
			}					
		}


		// ---------------------------------------------------------------------
		/// <summary>
		/// Add two matrices together and return the result as a new matrix 
		/// 
		/// Usage:
		///   Matrix m3 = Matrix.add (m1, m2);
		/// </summary>
		/// <param name="x">First Matrix</param>
		/// <param name="y">Second Matrix</param>
		/// <returns>Ruturns the sum of matrices, x and y</returns>
		// ---------------------------------------------------------------------
		public static Matrix add (Matrix x, Matrix y) {
			if (sameDimensions (x, y)) {
				Matrix result = new Matrix(x.r, x.c);
				for (int i=0; i<x.r; i++)
					for (int j=0; j<x.c; j++)
						result[i,j] = x[i,j] + y[i,j];
				return result;
			} else
				throw new EMatrixException ("Matrices must be the same dimension to perform addition"); 
		}
	    
	    
		// ---------------------------------------------------------------------
		/// <summary>
		/// Add a matrix to the current matrix. The current matrix is overwritten. 
		/// 
		/// Usage:
		///    m.add (m1)
		/// </summary>
		/// <param name="x">Right hand side of sum</param>
		// ---------------------------------------------------------------------
		public void add (Matrix x) {
			if (sameDimensions (this, x)) {
				for (int i=0; i<x.r; i++)
					for (int j=0; j<x.c; j++)
						this[i,j] = this[i,j] + x[i,j];
			} else
				throw new EMatrixException ("Matrices must be the same dimension to perform addition"); 
	       
		}
	    
		// ---------------------------------------------------------------------
		/// <summary>
		/// Static method to add a scalar value to a matrix and return a new matrix.
		/// 
		/// Usage:
		///    Matrix m1 = Matrix.add (m2, Pi);
		/// </summary>
		/// <param name="x">Matrix to modify</param>
		/// <param name="k">Scalar value</param>
		/// <returns>Returns m1 + k</returns>
		// ---------------------------------------------------------------------
		public static Matrix add (Matrix x, double k) {
			Matrix result = new Matrix (x.r, x.c);
			for (int i=0; i<x.r; i++)
				for (int j=0; j<x.c; j++)
					result[i,j] = x[i,j] + k;
			return result;
		}


		// ---------------------------------------------------------------------
		/// <summary>
		/// Add a scalar value to a matrix. The current matrix is overwritten.
		/// 
		/// Usage:
		///   m.add (3.1415)
		/// </summary>
		/// <param name="k">Scalar value</param>
		// ---------------------------------------------------------------------
		public void add (double k) {
			for (int i=0; i<r; i++)
				for (int j=0; j<c; j++)
					this[i,j] = this[i,j] + k;
		}


		// ---------------------------------------------------------------------
		/// <summary>
		/// Static method to subtract two matrices and return the result as a new matrix. 
		/// 
		/// Usage:
		///   Matrix m3 = Matrix.sub (m1, m2);
		/// </summary>
		/// <param name="x"></param>
		/// <param name="y"></param>
		/// <returns>Returns x - y</returns>
		// ---------------------------------------------------------------------
		public static Matrix sub (Matrix x, Matrix y) {
			if (sameDimensions (x, y)) {
				Matrix result = new Matrix(x.r, x.c);
				for (int i=0; i<x.r; i++)
					for (int j=0; j<x.c; j++)
						result[i,j] = x[i,j] - y[i,j];
				return result;
			} else
				throw new EMatrixException ("Matrices must be the same dimension to perform addition"); 
		}
	    
	    
		// ---------------------------------------------------------------------
		/// <summary>
		/// Subtract a matrix from the current matrix. The current matrix is overwritten. 
		/// 
		/// Usage:
		///    m.sub (m1);
		/// </summary>
		/// <param name="x"></param>
		// ---------------------------------------------------------------------
		public void sub (Matrix x) {
			if (sameDimensions (this, x)) {
				for (int i=0; i<x.r; i++)
					for (int j=0; j<x.c; j++)
						this[i,j] = this[i,j] - x[i,j];
			} else
				throw new EMatrixException ("Matrices must be the same dimension to perform subtraction"); 
		}

	    
		// ---------------------------------------------------------------------
		/// <summary>
		/// Subtract a scalar value from a matrix and return a new matrix. 
		/// 
		/// Usage:
		///   Matrix m1 = Matrix.sub (m2, 2.718);
		/// </summary>
		/// <param name="x">Matrix</param>
		/// <param name="k">Scalar</param>
		/// <returns>Returns (m2 - scalar)</returns>
		// ---------------------------------------------------------------------
		public static Matrix sub (Matrix x, double k) {
			Matrix result = new Matrix (x.r, x.c);
			for (int i=0; i<x.r; i++)
				for (int j=0; j<x.c; j++)
					result[i,j] = x[i,j] - k;
			return result;
		}
        
    
		// ---------------------------------------------------------------------
		/// <summary>
		/// Subtract a scalar value from the current matrix. The current matrix is overwritten. 
		/// 
		/// Usage:
		///   m.sub (3.1415);
		/// </summary>
		/// <param name="k">Scalar value</param>
		// ---------------------------------------------------------------------
		public void sub (double k) {
			for (int i=0; i<r; i++)
				for (int j=0; j<c; j++)
					this[i,j] = this[i,j] - k;
		}


		// ---------------------------------------------------------------------
		/// <summary>
		/// Multiply the current matrix by a scalar. The current matrix is overwritten.
		/// 
		/// Usage:
		///   m.mult (3.1415);
		/// </summary>
		/// <param name="k">Scalar value</param>
		// ---------------------------------------------------------------------
		public void mult (double k) {
			for (int i=0; i<this.r; i++)
				for (int j=0; j<this.c; j++)
					this[i, j] = this[i,j] * k;
		}


		// ---------------------------------------------------------------------
		/// <summary>
		/// Multiply the current matrix by a matrix. The current matrix
		/// is overritten and may change dimensions.
		/// 
		/// Usage:
		///   m.mult (A);
		/// </summary>
		/// <param name="A">Matrix</param>
		// ---------------------------------------------------------------------
		public void mult (Matrix A) {

			Matrix result = Matrix.mult (this, A);
			this.resizeSelf (result.r, result.c);
			result.copyTo (this);
		}


		// ---------------------------------------------------------------------
		/// <summary>
		/// Static method to multiply two matrices together. 
		/// 
		/// Usage: 
		///   Matrix C = Matrix.mult (A, B);
		/// </summary>
		/// <param name="m1">Left matrix</param>
		/// <param name="m2">Right matrix</param>
		/// <returns>Return the new matrix</returns>
		// ---------------------------------------------------------------------
		public static Matrix mult (Matrix m1, Matrix m2) {
		
			Matrix result = new Matrix (m1.r, m2.c);
			if (m1.c == m2.r) {
				for (int i=0; i<m1.rows; i++)
					for (int j=0; j<m2.cols; j++) {
						for (int k=0; k<m1.cols; k++)
							result[i,j] += m1[i,k] * m2[k,j];
					}
				return result;
			}
			else
				throw new EMatrixSizeError ("Incompatible matrix operands to multiply");
		}

		// ---------------------------------------------------------------------
		/// <summary>
		/// Computes the transpose of the current matrix. The current matrix is overwritten. 
		/// 
		/// Usage: 
		///   m.transposeSelf();
		/// </summary>
		// ---------------------------------------------------------------------
		public void transposeSelf () {

			int ir = r, ic = c;
			Matrix m = new Matrix (ic, ir);
			for (int i=0; i<ir; i++)
				for (int j=0; j<ic; j++)
					m[j,i] = this[i,j];
			this.mxx = m.mxx;
			this.r = ic;
			this.c = ir;
			
			// Swap the row and column labels
			string[] tmp = new string[ir];
			rowNames.CopyTo (tmp, 0);
			rowNames = new string[ic];
			columnNames.CopyTo (rowNames, 0);
			columnNames = new string[ir];
			tmp.CopyTo (columnNames, 0);
		}
	    

		// ---------------------------------------------------------------------
		/// <summary>
		/// Computes the transpose of the current matrix and returns it
		/// The current matrix is not changed
		/// 
		/// Usage:
		///   A = Matrix.transpose(B);
		/// </summary>
		// ---------------------------------------------------------------------
		public static Matrix transpose (Matrix A) {

			int ir = A.r; int ic = A.c;
			Matrix result = new Matrix (ic, ir);
			for (int i=0; i<ir; i++)
				for (int j=0; j<ic; j++)
					result[j,i] = A[i,j];

			result.rowNames = new string[A.c];
			result.columnNames = new string[A.r];

			A.rowNames.CopyTo (result.columnNames, 0);
			A.columnNames.CopyTo (result.rowNames, 0);
			return result;
		}

		
		// ********************************************************************
		/// <summary>
		/// Determines if a matrix is non-singular or not. Returns true if it is 
		/// non-singular. Singularity is determined by checking whether the 
		/// rank (m) is less than the number of rows. 
		///  
		/// Usage:
		/// <code>
		/// if (m.isNonSingular()) ....
		///  </code>
		/// </summary>
		/// <returns>True if the matrix is non-singular</returns>
		// ********************************************************************
		public bool isNonSingular() {
			int n = this.c;
			if (this.rank() < this.r)
				return false;
			else return true;
		}

		
		// ********************************************************************
		/// <summary>
		/// Solve finds and returns the solution x to the linear system, A*x = B, 
		/// using LU decomposition. B is usually a column vector, yielding a column 
		/// vector solution, x. However, as in many implementations, B may have 
		/// multiple columns such that, each column of B generates a solution vector 
		/// in the solution matrix, x. The inverse of a matrix A can 
		/// therefore generated by setting B to the identity matrix.
		///  
		/// Usage:
		/// <code>
		///   solution = A.solve (B)
		///  </code>
		/// </summary>
		/// <returns>Returns the solution vectors</returns>
		// ********************************************************************
		public Matrix solve (Matrix B) {
			
			int m = this.r;
			int n = this.c;

			if (B.r != m) {
				throw new System.ArgumentException("Matrix row dimensions must agree.");
			}

			// Compute the LU decomposition
			LUDecomposition LUD = new LUDecomposition (this);

			// Check if the matrix is singular
			for (int j = 0; j < n; j++) {
				if (LUD.LU[j,j] == 0)
					throw new System.SystemException("Matrix is singular.");
			}
			
			// Copy right hand side with pivoting
			int nx = B.c;
			Matrix X = new Matrix(LUD.pivotVector.Length, (nx-1) + 1);
			for (int i = 0; i < LUD.pivotVector.Length; i++) {
				for (int j = 0; j <= nx-1; j++) {
					X[i,j] = B[LUD.pivotVector[i],j];
				}
			}

			// Solve L*Y = B(piv,:)
			for (int k = 0; k < n; k++) {
				for (int i = k + 1; i < n; i++) {
					for (int j = 0; j < nx; j++) {
						X[i,j] -= X[k,j] * LUD.LU[i,k];
					}
				}
			}
			// Solve U*X = Y;
			for (int k = n - 1; k >= 0; k--) {
				for (int j = 0; j < nx; j++) {
					X[k,j] /= LUD.LU[k,k];
				}
				for (int i = 0; i < k; i++) {
					for (int j = 0; j < nx; j++) {
						X[i,j] -= X[k,j] * LUD.LU[i,k];
					}
				}
			}
			return X;
		}

		
		// ********************************************************************
		/// <summary>
		/// Compute the inverse of a matrix. Inverse is returned by the call.
		/// The current matrix is not modified.
		///    
		/// Usage:  
		///   inv = A.invert()     
		/// </summary>
		/// <returns>Returns the inverse</returns>
		// ******************************************************************** 
		public Matrix invert() {

			if (this.r != this.c)
				throw new EMatrixSizeError ("Matrix must be square in order to compute the inverse");

			return this.solve (Matrix.identity(this.r));
		}


		// ********************************************************************
		/// <summary>
		/// Compute the inverse of the current matrix                            
		/// The current matrix is overwritten with the inverse  
		///    
		/// Usage:  
		///   A.invertSelf()  
		/// </summary>
		// ******************************************************************** 
		public void invertSelf() {
			Matrix cpy = this.invert();
			cpy.copyTo (this);
		}


		// ******************************************************************** 
		/// <summary>
		/// Find the determinant of a square matrix
		/// </summary>
		/// <param name="A">Find the determinant of the matrix</param>
		/// <returns>Returns a scalar value equal to the determinant</returns>
		// ******************************************************************** 
		public static double det (Matrix A) {

			if (A.r == A.c) {
			    LUDecomposition LUD = new LUDecomposition (A);
				double result = LUD.pivotSign;
				for (int i=0; i<A.r; i++)
					result = result * LUD.LU[i,i];
				return result;
			} else
				throw new EMatrixSizeError ("Determinant: Matrix must be square");
		}


		// ******************************************************************** 
		/// <summary>
		/// Find the determinant of the current square matrix
		/// </summary>
		/// <returns>Returns a scalar value equal to the determinant</returns>
		// ******************************************************************** 
		public double det () {

			if (this.r == this.c) {
				LUDecomposition LUD = new LUDecomposition (this);
				double result = LUD.pivotSign;
				for (int i=0; i<this.r; i++)
					result = result * LUD.LU[i,i];
				return result;
			} else
				throw new EMatrixSizeError ("Determinant: Matrix must be square");
		}

		// ******************************************************************** 
		/// <summary>
		/// Returns the rank of the current matrix. Uses rowEchelon to compute the value
		/// </summary>
		/// <returns>Rank</returns>
		// ******************************************************************** 
		public int rank() {
			int rank = 0;
			this.rowEchelon (out rank);
			return rank;
		}


		// ******************************************************************** 
		/// <summary>
		/// Reduce the current matrix to the fully reduced row echelon form. 
		/// The current matrix is overwritten in the process
		/// From James D. Reilly, Fort Wayne, Indiana
		/// </summary>
		/// <param name="rank">Rank of the matrix</param>
		/// <returns>Returns the determinant</returns>
		// ******************************************************************** 
		public double rowEchelon(out int rank) {
            int x;
            double Pivot, Temp;
            double Determinant;
            int PivotRow, PivotCol;

            PivotRow = 0;
            PivotCol = 0;
            Determinant = 1.0;
            do {
                // Find largest pivot.  Search for a number below
                // the current pivot with a greater absolute value.
                x = PivotRow;
                for(int C=PivotRow; C<rows; C++)
                    if(Math.Abs(this[C, PivotCol]) > Math.Abs(this[x, PivotCol]))
                        x = C;

                if(x != PivotRow) {
                    // If here, there is a better pivot choice somewhere
                    // below the current pivot.
                    // Interchange the pivot row with the row
                    // containing the largest pivot, x
                    for(int b=0; b<cols; b++) {
                        Temp = this[PivotRow, b];
                        this[PivotRow, b] = this[x, b];
                        this[x, b] = Temp;
                    }
                    Determinant = -Determinant;
                }

                Pivot = this[PivotRow, PivotCol];
                Determinant = Pivot*Determinant;

                if(Pivot != 0.0) {
                    // Introduce a '1' at the pivot point
                    for(int b=0; b < cols; b++)
                        this[PivotRow, b] = this[PivotRow, b]/Pivot;
						     
                    for(int b=0; b < rows; b++) {
                        // Eliminate (make zero) all elements above
                        // and below the pivot.
                        // Skip over the pivot row when we come to it.
                        if(b != PivotRow) {
                            Pivot = this[b, PivotCol];
                            for(int C = PivotRow; C < cols; C++)
                                this[b, C] = this[b, C] - this[PivotRow, C]*Pivot;
                        }
                    }
                    PivotRow++; // Next row
                }
                PivotCol++;  // Next column
            }
            while((PivotRow < rows) && (PivotCol < cols)); // Reached an edge yet?

            // Finally compute the rank
			double[] sq = new double[this.r];
			for (int i = 0; i < this.r; i++)
				for (int j = 0; j < this.c; j++)
					sq[i] = sq[i] + Math.Abs (this[i,j]);
			rank = 0;
			for (int i = 0; i < this.r; i++)
                if (sq[i] < 1E-6)
					rank++;
			return Determinant;
        }

	}
}
